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FOREWORD 

This is the final report of a theoretical investigation of 
M S-IC Base Radiation and Plume Interaction Reversed Flow 
Study. " This report presents documentation for the base flow 
analysis effort performed on the study. Documentation for 
other efforts have been published previously and are not repeated 
here (NASA CR-61321, "User's Manual for "RAVFAC" - A 
Radiation View Factor Digital Computer Program," by J. K. 

Lovin and A. W. Lubkowitz, LMSC/HREC D 148620, November 
1969). This work was performed by Lockheed Missiles & Space 
Company, Huntsville Research & Engineering Center, for the 
Aero-Astrodynamics Laboratory of the National Aeronautics and 
Space Administration, George C. Marshall Space Flight Center, 
under Contract NAS8-30154. 

Technical monitors for this contract were Mr. H. B. Wilson, 
Jr., Mr. R. F. Elkin, and Mr. E. B. Brewer of the Thermal 
Environment Branch, Aero-Astrodynamics Laboratory, MSFC. 
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SUMMARY 

Problems associated with computing the three-dimensional base re- 
circulation flow of a perfect gas from multiple-clustered nozzles are dis- 
cussed. A Lax-Wendroff time -dependent, finite-difference numerical scheme 
is presented as the most desirable technique for computing this flow field. 
Results of a preliminary investigation indicated that including viscous terms 
in this analysis was impractical because of the number of grid points required 
for accurate resolution within strong viscous regions. However, inclusion 
of viscous terms was deemed unnecessary for computing the overall flow field. 
A cylindrical coordinate system with the axis at the center of the nozzle cluster 
simplified the specification of some boundary conditions, but extrapolation was 
required to specify conditions downstream of the nozzles and around the exterio 
of the nozzle wall. Stretching transformations are presented which concentrate 
grid points in the regions of the plume interaction. A possible modification to 
the numerical scheme which would permit using a steady form of the energy 
equation and, hence, reduce computer storage requirements is discussed. 
Stability of the numerical technique used in this study is also analyzed. Two 
techniques are presented for controlling the time step to ensure stability during 
the computation. The first technique monitors the total change between time 
steps of all dependent variables over the entire flow field and controls the time 
step to ensure that this total change decreases. The second technique prevents 
using a time step large enough to cause a provisional value of density to go 
negative at any grid point. Conclusions resulting from this analysis and recom 
mendations for future investigations are presented. 
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NOMENCLATURE 

a centrifugal force vector for cylindrical coordinate system (Eq. (2*10e) ) 

b a r/ae 

c speed of sound 

c v specific heat at constant volume 

E total energy 

e internal energy per unit mass 

f vector used to compute changes in x direction (Eq, (2. 10b) ) 

g vector used to compute changes in r direction (Eq. (2,10c) ) 

g quantity used to monitor change of dependent variables (Eq. (2.28) 

h vector used to compute changes in 0 direction (Eq, (2, lOd) ) 

stretching constant in r direction 
K t stretching constant in T direction 

n time step 

NCN number of concentric nozzles 

p pressure 

q velocity vector 

R axisymmetric radial distance 

& gas constant 

r radial distance (Fig. 2) 

r" incremental parameter for stretching in r direction 

T temperature 
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NOMENCLATURE (Continued) 


t 

U 

u 

V 

V 
w 
X 

X 

X 

o 

Greek 

P 

y 

A 

6 

e 

e 

o 

p 

T 

? 

4 

u 


time 

axisymmetric velocity component in X direction (Fig. 4) 
velocity component in x direction 

axisymmetric velocity component in R direction (Fig. 4) 
velocity component in r direction 
velocity component in 0 direction 

axisymmetric longitudinal distance aft of nozzle exit plane 
longitudinal distance aft of nozzle exit plane 

negative longitudinal distance from nozzle exit plane to heat shield 
radial distance from nozzle centerline to centerline of cluster 


tangency angle of nozzle exterior wall in r-0 plane (Fig. 3) 

ratio of specific heats 

step size 

spatial step size 

angular direction 

angle defining reflection plane between nozzles (Fig. 2) 
density 

normalized parameter in 0 direction 

incremental parameter for stretching in 0 direction 

angle between axisymmetric velocity vector and nozzle centerline 
(Fig. 4) 

vector form of dependent variables (Eq. (2.10a) ) 
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NOMENCLATURE (Continued) 


Subscripts 

w condition at exterior of nozzle wall 

t total reference condition 

i,j,k indices for grid points in x, r, 0 directions, respectively 

0 value at previous time 

1 value at current time 


vii 
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Section 1 
INTRODUCTION 

Interaction of exhaust plumes from multiple rocket nozzle configurations, 
as shown in the example configuration of Fig. 1, creates a complex three- 
dimensional network of shock waves. If the impingement angle is large, 
which is usually the case of highly underexpanded plumes, these shock waves 
may detach from a plane of symmetry existing between the plumes. This 
will allow extremely hot gases to recirculate around the heat shield and the 
exterior of the rocket nozzle; thus, heat transfer to these components could 
be significantly increased. Although base recirculation is a well-recognized 
problem of launch and space vehicles, the three-dimensionality of the flow 
has prevented thorough analyses. Experimental investigations (Ref. 1) have 
yielded some insight into the mechanism of these flows, but applicability of 
these data are very limited. The purpose of this study is to develop a theo- 
retical technique which would provide the required flexibility for investigating 
the base recirculation flow field of multi-nozzle configurations. 
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Fig. 1 — Example of Four Nozzle Configurations 
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Section 2 

TECHN IC A L D ISC USSION 

2.) SELECTION OF TECHNIQUE FOR SOLUTION OF FLOW FIELD 
EQUATIONS 

Of all the technique used to solve fluid mechanic problems, cloted- 
torm analytical solutions are obviously the rnoet desirable. Unfortunately, 
the conservation equations of fluid motion contain nonlinear partial differen- 
tial equations. Closed-form solutions to these equations require that they 
be linearized either by transformations or by simplifying assumptions. The 
co^npiexity of the equations describing the three-dimensional base recircu- 
lation flow field, which contains shock waves and resulting regions of sub- 
sonic and supersonic flow, precludes either transformations or approximations 
that are valid over the entire field. 

Numerical solutions of the conservation equations of this recirculation 
flow present the only alternative. Spatial marching techniq”es, such as the 
method of characteristics or the steady-state form of the Lax-Wendroff 
technique, have been used with apparent success for plume interaction prob- 
lems which are limited to oblique shock waves at the interacting reflection 
planes (Refs. 2 and 3). Such initial-value techniques are restricted to 
regions of hyperbolic flow; hence, they are entirely inadequate for analysis 
of flow containing detached shocks with regions in which the equations are 
elliptic i . nature. 

Two classes of numerical techniques which appear most attractive for 
the solution of this type of problem are asymptotic time-dependent unsteady 
techniques and relaxation techniques. Both classes involve solving finite- 
uifference approximations to the Eulerian equations of motion. The time- 
dependent techniques use the principle that the unsteady Eulerian equations 
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are hyperbolic in nature; thus, marching forward in time is possible. If 
the time step is appropriately small such that the numerical scheme is 
stable, the unsteady flow will asymptotically approach the steady-state 
solution. 

One of the few formalized relaxation techniques for fluid mechanics 
applications was recently presented by Prozan (Ref. 4). His technique, 
known as the Error Minimization Technique (EMT), uses the steady Eulerian 
form of the conservation equations, but these equations are set equal to an 
error term for each equation at each grid point. Because an initial flow field 
with imposed boundary conditions does not satisfy these conservation equa- 
tions, this error term is non-zero. A positive definite merit function of the 
error term is defined, and the merit function is decreased through descent 
techniques. As this merit function asymptotically approaches zero, the flow 
variables approach their correct steady-state values. Comparisons with 
time -dependent techniques indicate that the EMT is equally as accurate and 
requires comparable computation time as the unsteady techniques require 
for convergence. 

The feasibility of including viscous terms in an analysis of this type 
was considered. Both the unsteady and relaxation methods are capable of 
solutions involving these terms. However, a preliminary investigation indi- 
cates that a line mesh grid system was required to obtain resolution in re- 
gions where the viscous stresses dominate. When considering the staggering 
number of grid points required to define the inviscid three-dimensional flow 
field, it was concluded that the inclusion of enough additional grid points to 
yield meaningful results of a viscous flow field of this type was impractical. 

It was further concluded that the exclusion of viscous forces from the analysis 
would not significantly affect the characteristics of the overall flow field, 
although the possibility of rationalization is conceded. 

Infinitesimally thin shock waves of various strength will obviously exist 
throughout the flow field under consideration. State-of-the-art, finite- 
difference techniques are not able to correctly isolate the location of shocks 
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between grid points, treat them as discrete discontinuities, and apply shock 
jun.p conditions. An attempt to improve the state of the art to accomplish 
this was unsuccessful. An alternative approach was to introduce artificial 
viscosity into the equations. This may be done explicitly as presented by 
von Neumann and Richtmyer (Ref. 5) or implicitly as done by Lax and Wendroff 
(Ref. 6). Either approach introduces dissipation terms which become large 
in regions of high compression. As a result, shock waves occur over several 
grid points. The effect of this smearing can often be reduced by using a finer 
grid mesh. 

Presumably, because of the similarity of the formulation of the equa- 
tions, shock smearing versions could exist for either EMT or unsteady tech- 
niques. However, several versions of the unsteady techniques have proven 
quite capable of computing flow fields with shock waves; whereas a shock 
smearing version of EMT would require additional development. For this 
reason, a i ecent version of the finite -difference, two-step Lax-Wendroff 
unsteacj technique is presented in Section 2.4 of this report for solving the 
three-dimensional base recirculation flow field. This technique has been 
coded in FORTRAN V for use on the Univac 1108 digital computer. The 
following sections of this report describe the application of this technique 
to the base recirculation problem, as well as a derivation of the finite- 
difference scheme. 

2. 2 GEOMETRY CONSIDERATIONS 

The computation technique has much versatility in respect to geometry 
considerations. Flow from any number of nozzles may be computed so long 
as they are positioned concentrically about a central point. Figure 2 pre- 
sents an example of a four-nozzle configuration. The end view reveals eight 
reflection planes for this configuration. Because of this symmetry, it is 
necessary to investigate only one -eighth of the entire flow field for this con- 
figuration. All distances are nondimensionalized by the exit diameter of the 
nozzle. The effect of nozzle separation distance relative to the exit diameter 
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may be studied by choosing various values of y . The heat shield location 

is designated x , where x <0. 

° o o — 

2.2.1 Coordinate System 

The choice of a cylindrical coordinate system is natural when con- 
sidering the region of analysis. The angle from the reference line to the 
first plane of symmetry 0 q is computed from the geometry by 


JT 7T 

o = 2 ‘ NCN 


( 2 . 1 ) 


where NCN is the number of concentric nozzles. For convenience of nota- 
tion, a normalized parameter T is defined as 

0 

- 2 - ; (2 
0 

o 

thus, T = 0 and 1 represents the two reflection planes. 

2.2.2 Exterior Nozzle Wall Equations 

The exterior wall of the nozzle is represented as a cylinder as shown 
in Fig. 2. For this offset cylindrical coordinate system, the equation for 
this nozzle exterior wall is 



sin0 

w 


2^2 1 
r + y - -r 
w 7 o 4 

2y r 
7 o w 


(2.3) 


if (r - |) < y Q < (r + j) and x < 0 . 

By defining 3 w as the angle between the tangent of the exterior wall 
and the horizontal in the r - 0 plane (Fig. 3), then 
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tan P w 


dr 

Te 


w 


w 


= y cos0 
7 o w 


y 


o 


sinG 

w 



2 

y o 


2 

COS 



(2.4) 


This angle will be used to specify tangential flow about the nozzle exterior 
wall. 


2.2.3 Stretching Transformations 

By visualizing flow for the nozzles in Fig. 2, it may be seen that plumes 
will first interact and cause shock waves along T = 0 . Slightly farther down- 
stream all four plumes will intersect at r = 0 causing a complex flow 
pattern to occur inside a diamond-shape shock wave for the four-nozzle 
configuration. For these reasons, increased emphasis should be placed 
on defining the region near T = 0 and r = 0. The following stretching trans- 
formation is suggested to provide a fine mesh grid near T = 0 : 

T = T - K r sin(7rr) (2.5) 

where 

K r = arbitrary constant less than — 

T = parameter varied from 0 to 1 by even increments 
A similar stretching transformation in the r direction is suggested as 

r = tan(K r 'r) (2.6) 

where 

K r = arbitrary constant 

r* = parameter which may be varied from 0 to 7 r / 2K f 
by even increments 
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2.2.4 Number of Grid Points 


Approximately 50 grid points in the x direction, 25 in the r direction, 
and 2 3 in the T direction are estimated to be the minimum number of grid 
points required for adequate resolution of a typical problem for the geometry 
shown in Fig. 2. This total comes to nearly 30,000 grid points. It will be 
shown in Section 2.4 that 10 dependent variables must be stored at each grid 
point. This requires storing 300,000 variables which greatly exceeds the 
core storage of any digital computer available locally. A rapid access drum 
storage routine of the Univac 1108 digital computer provides an alternative 
to core storage. This routine known as NTRAN has been modified to re- 
linquish control of the machine while performing the time-consuming chore 
of repositioning the drum to a desired location. Combining a restart capability 
with this routine would be desirable to prevent loss of converging flowfield 
parameters through machine error. 

2. 3 BASIC DIFFERENTIAL EQUATIONS 


The unsteady conservation of mass and energy equations, the three un- 
steady Euler equations, and the perfect gas relation give six equations in six 
unknowns. These unknowns are the three velocity vectors (u, v, w), pressure 
(p), total energy (E), and density (p). The pressure may immediately be 
eliminated from these equations by considering the perfect gas relation 


A 

P 


A A A 

pSt T 


AAA 

= (y - 1) pC v T 


= (y - 1) 


aa 

P e 


Nondimensionalizing by the reference total density p^ and speed of sound 

A 

c^ gives 


£ 

Note: Quantities with ( A ) denote dimensional variables. 
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P = (y - 1) pe 


Using the definition of the total energy as the sum of the internal and kinetic 
energies 


E .d['*T’ 2 ] 


(2.7) 


where 


q 2 * 


2 2 
u + v 


+ w 


yield the pressure as 

P = (y - 1) |e - ^J- , (2.8) 

This expression is used to eliminate pressure from the Euler equations and 
the conservation of energy equation. This leaves the following set of five 
equations and five unknowns shown in vector notation as 


where 


9ui 

w = 


f8f 


, J_ 8h 
r r 


f] 


(’\ 

pv 

V>7 


(2.9) 


(2.10a) 
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(2.10b) 



(pw) 

(pu) \pw ) 

(pv) (pw) 

(v . 1)E . (i) (pw) 2 - (^) (|) [(pu) 2 

y (»w) - [4-4 (-y) <pw) [(pu ) 2 + (pv ) 2 1 

p 




(2. lOd) 
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(pv) 


p ) (pv) { 


yE 


(pu) (pv) 

(y) [(P v > 2 - <P W > 2 ] 

(pv) (pw) 

(^T^) (t) [ (pu)2 + (pv) ' + (pw)2 ]} 


(2. lOe) 


These equations, expressed in conservative form, treat (pu), (pv), and (pw) 
as dependent variables instead of u, v, w. 


An alternate form of these equations in which E may be expressed as 
a function of the remaining four dependent variables is proposed. Combining 
the adiabatic, perfect gas relation, 



with the specific internal energy relation. 


A 

e 


A 

c 


V 


A 

T 


A 

ft 

" (y - l) 


A 

T 


and the thermally perfect expression of the speed of sound, 


A2 

c* 


y«t T t 
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gives 


= i(- 

V \i 


1 

TT 



The total energy E may then be explicitly expressed as a function of the other 
four dependent variables by 



( 2 . 11 ) 


Equation (2. 11), which is not in diffe?rential form, can be used to decrease the 
computer storage variables by 20%. Preliminary investigations on simplified 
flow fields indicate that convergence time is not significantly altered by using 
Eq. (2.11) to solve for the energy directly. The computer program does not 
include this option because of the preliminary nature of these results. 

2.4 SOLUTION OF UNSTEADY FINITE-DIFFERENCE EQUATIONS 


The solution to Eq. (2.9) begins by expanding the vector Eq. (2. 10a) in 
the Taylor series 


u 


1 

i. j» k 




At + C(At) 2 


( 2 . 12 ) 


where superscript t refer to time and subscripts i, j, k refer to grid point 
numbers in the x, r, 0 directions Direct substitution of Eq, (2.9) into Eq. (2.12) 
gives 
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Any finite-difference analogy to Eq. (2.13) is unconditionally unstable; however, 

this equation can be modified into a form which is conditionally stable. One 

popular technique is to retain higher order ttrms, but this involves the com- 

2 2 

plicated and computer time-consuming process of expressing 8 u/9t as 
a function of 9^f/9x^, 9^f/9x 9r , etc. 

A preferred method of stabilizing Eq. (2.13) is by using a two-step process. 
This process, known as the two-step Lax- Wendroff technique, computes 
provisional values of u in the first step which are then averaged with previous 
value* in a second stable step. A three-dimensional, cylindrical coordinate 
adaptation of a two-step version recently presented by MacCormack (Ref- 7) 
is suggested for the solution of the base recirculation flow equations. This 
technique was chosen because it ie easy to compute and has second-order 
accuracy. 


Consider the one -dimensional form of Eq. (2. 13„ 



i 



(H) At + 0(At)2 - 

' 'i 


(2.14) 


The central difference analogy of the partial derivative of Eq. (2. 14) may be 
expressed as 


/ 9 f \ ^ f i+l ‘ 

\JZ). zZT 


- f. 
1 


f. - 

l 


f i-i 


2 Ax 


+ 0(Ax) 


(2.15) 


The first step of the technique which computes provisional values, steps 
forward in time 2At and uses the forward portion of Eq. (2.15) or 


*1 l 0 IfO At 

u i = u i ' \ f i+i ' f i / a7 * 


(2.16) 


The second and stabilizing step of the process provides second-order accuracy 
by using the backward portion of E '. (2.15). This step marches forward in time 
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At from tg and averages provisional and original values of u by 


‘o . . .h 


At 
2 Ax 


*i u i + u i is s \ 

- 2 l £ i - h-i) 


- - ( £ i+i - £ ‘° + £ !‘ - £ *! i ) + °^*» 2 + < 2 - 17 > 


A three-d. nensional version of Eqs. (2.16) and (2.17) is presented 


First Stet 


n+l _ n At \(m _ ji ( n n \ 

U i,j,k U i,j,k " 6 | \ l+l, j* k ' i,j,k y g i, j+1, k “ g i, j, k/ 


Second Step: 


n+l _1_ ( n n+l At [/.n+l ,i +1 \ / n+l n+l \ 

J i,j,k ~ 2 ( i, j,k + U i,j,k 6 [\ i,j,k " i-l,j,k/ \ g i,j,k " g i,j-l,k/ 

+ 2. (h 577 - h* 77 .)] - ¥ (a 577 )! (2.1 

r \ i,j,k i,j,k-l/J r \ l.j.k/j 

vhere 

n = superscript used to increment t by t = (n+l) At 

i, j,k = subscripts referring to grid points in the x,r,r directions, 

respectively 

6 = Ax = Ar = AT 

9 T NCN 
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As suggested by MacCormack, increased stability (larger permissible 

values of At) might be obtained by a cyclic procedure using the 

forward and backward portions of Eq. (2.15) in the first and second step. 

Following this suggestion, Eq. (2.13), which consists of three separate 

3 

partial derivatives, was approximated by eight (i.e., 2 ) combinations of 
forward and backward finite difference equations. Equations (2.18) and (2.19) 
represent the first of these eight combinations. Programming these eight 
schemes was simplified by using a variable indices technique suggested by 
E. B. Brewer, NASA Technical Monitor, 

2.5 BOUNDARY CONDITIONS 

The choice of a cylindrical coordinate system (Fig, 2) simplified the 
specification of boundary conditions. Along the planes T = 0, 1 , reflection 
principles can be applied with no qualification. Axial symmetry provides 
the boundary condition at r = 0. Because of the singularity of the equations 
at r = 0 , the grid was chosen such that this line was straddled. Downstream 
boundary conditions for maximum values of x were obtained by linear extrapo- 
lation. Flow in this region is supersonic; therefore, this approximation appears 
warranted. Similarly, conditions at maximum values of r were obtained by 
linear extrapolation. These grid points of r which lie within the plv.me are 
again in a supersonic flow regime. It is anticipated that others (i.e., x < 0) 
will have a negligible effect on the plume interaction region. 

Conditions along the heat shield are also easy to specify. As pointed out 
by Moretti (Ref. 8), applying reflection techniques at solid boundaries might 
constitute improperly forcing to zero certain derivatives in the direction 
normal to the wall. This could bias the solution. However, it would be 
possible to substitute a mirror image flow neld behind the flat surface of the 
heat shield without altering the flow. Thus, reflection principles along the 
heat shield are applicable and were used. 

Specifying boundary conditions to ensure that flow about the exterior 
wall of the nozzle is tangential is more difficult. Figure 3, which shows 
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the exterior of the nozzle in the r - 0 plane, reveals that grid points will 
not lie on the surface of this wall. From geometry considerations specified 
in Eq. (2.3), three simple tests will reveal if a grid point lies inside of the 
nozzle plane. Interior grid points adjacent to the wall may then be used to 
specify boundary conditions at the wall. An extrapolation in the T (i.e., 0) 
direction using values of p, pu, and E at neighboring exterior points pro- 
vide values of p, pu, and E at the wall and interior adjacent points. Re- 
arranging Eq. (2.11) to solve for q gives a wall value of 

2 _ 2 V | Ew 1 1 r, , 

^ Y " 1 L P w " Y<V - DJ * 

Multiplying the definition 

2 _ 2 2 2 
q = u + v + w 


by p vields 
1 w 


? 7 2 2 2 

(P q. ) - (pu) C + p (v + w^ ) 
'^w ^w w w w w 


or which may be rearranged as 


2,2 2 , 

p (v + w ) 
w w w 


« > w' i w |2 - <pu) w 


Applying the wall tangency condition gives 


<PV) W = - ( ' ,u) t] s “ |3 w 


and 


(P w). 


w 




where may be computed from Eq. (2.4) 


( 2 . 21 ) 


( 2 . 22 ) 
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Boundary values of pv and pw at adjacent interior grid points are then 
obtained by extrapolating in the T direction using values of pv and pw at 
exterior grid points and those computed at the wall by the tangency condition. 
Admittedly, this technique involves approximations and assumptions, but it 
is felt that more sophistication is unwarranted when considering computation 
time. 


Conditions at the nozzle exit are input as initial conditions and not 
allowed to vary. Grid points at this exit may also be located by using Eq. (2.3). 

2.6 INITIAL CONDITIONS 

It is desirable to input as initial conditions the best possible values 
of the dependent variables at each grid point throughout the field. This 
should significantly decrease the computation time required for convergence 
to the steady-state solution. For this reason, provisions were made in the 
computer program for inputting results of a method-of-characteristics (MOC) 
solution of a single-nozzle configuration. 

An interpolation routine is used to determine flowfield parameters at 
grid points which will not coincide with the axisymmetric characteristic grid 
points. Figure 4 illustrates the relation between the geometry of the axi- 
symmetric flow of the nozzle and the coordinate system used in this analysis. 
The calculation of p and pu are obtained directly from the interpolation 
routine which gives q, p, and at any specified point. The energy E may 
then be computed from Eq. (2.11). 

The two remaining flow variables, pv and pw, must be computed from 
geometry considerations. At a given grid point, values of r and 0 are 
known, and the axisymmetric radial distance may be computed by 

R = ^r^ - 2ry Q sin0 + y^ . (2.23) 
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The axisymmetric radial velocity at that point is given by 


V = q sin<£ . 


(2.24) 


Using Eqs. (2.23) and (2.24) and the fact that the flow is axisymmetric give 


and 


Pv 


(r - y Q sinG) pV 
R 


(2.25) 


p w 


- ~ pv cose 


(2.26) 


For grid points lying outside of the axisymmetric plume, velocity 
components are set equal to zero and ambient conditions are prescribed 
for p and E . 


Boundary conditions along the T = 0 plane and the r = 0 line are the 
only ones not identically satisfied. The reflection plane T = 0 is the plane 
of impingement of the two plumes; hence, care must be taken in specifying 
initial conditions along this plane which will not cause an initial instability. 
Following a suggestion by D'Attorre (Ref. 9), density and energy will be left 
as computed by the MOC program. Likewise, the computed parameter pu 
will be left unchanged, but pw will be set to zero with the calculated value 
of pw added to the calculated value of pu at that point. 

Reflection principles are used to specify boundary values of parameters 
across the T= 0 plane, and sxial symmetry is used across the r = 0 line. 

2.7 STABILITY 

A stable time -dependent technique is one which begins with an initial 
unsteady flow field and modifies flow variables so that they asymptotically 
converge to steady-state values. The maximum possible time step is given 
by 
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(2.27) 


which is the well-known Courant-Friedricks-Lewy (CFL,) condition. For 
this analysis, it should be noted that the appropriate 6 to be considered for 
Eq. (2.27) is the smallest value of Ax, ir, or A0 , regardless of stretching 
transformations used. 


A more restrictive stability condition for unsteady techniques is pre- 
sented by Von Neumann. Unfortunately, this condition can only be applied to 
linearized equations, b it linearization of the flow equations over the entire 
flow field is impos jibh*. because of the presence of shock waves and boundaries. 
It was, therefore, concluded that any attempt to predict a priori a stable time 
step, and claim that this time step should apply over the entire field through- 
out the analysis, would be unrealistic. Alternative approaches for measuring 
stability were, therefore, sought. 

One alternative approach is to use a technique similar to that used by 
Prozan in the EMT. By defining 



it is possible to monitor the total change of flow variables from one time 
step to the next. If it is appropriate to assume that the initial flow field, 
with boundary conditions applied, is more different from the steady-state 
flow than at any time during the convergence process, then g should be 
largest initially. By monitoring g and controlling the time step, it is 
possible to force g to monotonically decrease toward the asymptotic value 
of zero. Comparisons of g to a previous value should be made only after 
a complete cycle of the difference schemes discussed in Section 2.4 because 
some combinations of forward and backward differences are more stable 
and allow less change than others. It is possible that this technique of 
controlling the time step by monitoring g could result in neutral stability 
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(i.e. , instability prevented, but the solution does not converge). Results 
obtained by this technique should be analyzed to ensure that the solution 
represents a reasonable steady-state flow field. 

A second approach to monitoring stability is the antimatter (negative 
density) method. An initial time step of slightly less than that given by the 
CFL condition (Eq. 2.27) is assumed. Throughout the computer calculations, 
the provisional value of density calculated at each grid point during the first 
step is monitored. If this density should become negative, the time step is 
decreased and new provisional values are computed over the entire flow field. 
The basis for this procedure is the fact that all variables change more during 
the first step than during the second; hence, stability of variables computed by 
the second step might be maintained by rejecting an impossible value (negative 
density) computed during the first step. The danger of this approach is that 
the solution could have diverged to such an unrealistic condition when negative 
density occurs that recovery is very difficult. 


2-21 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 


LMSC/HREC D162322 


Section 3 
CONCLUSIONS 


The following conclusions were reached as a result of this analysis: 

1. Closed form analytical methods and spatial marching numerical 
techniques are inadequate for treating base recirculation flows 
with regions of subsonic flow; however, either time-dependent 
or relaxation finite difference numerical techniques may be used 
to analyze these flows. 

2. Techniques are not currently available for locating and treating 
shock waves as discrete discontinuities between finite difference 
grid points; therefore, shock smearing techniques using pseudo 
viscosity must be used. 

3. It appeared ..either practical nor necessary to include viscous 
terms in this analysis. 

4. A recent version of a two-step Lax-Wendroff time -dependent, 
finite -difference technique appeared to be the most attractive 
numerical scheme for solving the three-dimensional base 
recirculation flow field. 

5. A cylindrical coordinate system greatly simplified the appli- 
cation of boundary conditions. 

6. Stretching transformations are desirable to concentrate grid 
points near regions where the flow is complicated by shock 
waves of plume interactions. 

7. A rapid access drum storage routine of the Univac 1108 digital 
computer provides an acceptable alternative to the limiting 
core storage of this computer. 

8. It is possible to solve for the total energy of the flow at each 
grid point as a function of the local density and velocity com- 
ponents. This constitutes using a steady- state form of the 
energy equation and permits reduction of storage variables 
by 20%. 

9. Computation time can be significantly reduced by beginning 
with an initial flow field from a single nozzle calculated by 
the method of characteristics. 
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10. Because of the complexity of the flow field, it is not possible 
to mathematically calculate the largest possible time step 
which will yield a stable solution. An alternative technique 
is suggested by which the maximum stable time step is esti- 
mated as slightly less than that given by the CFL condition. 
Stability is then monitored either by comparing changes of 
flow variables or checking for negative density. 
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Section 4 

RECOMMENDATIONS 

Although the objective of this investigation was to develop a technique 
for predicting the base recirculation flow field of Saturn-class launch vehicles, 
the method appears applicable for Space Shuttle configurations. It is, there- 
fore, recommended that modifications to this technique be explored so that 
base recirculation of Space Shuttle vehicles may be computed. 

In principle, the numerical technique described in Section 2.4 of this 
report should be capable of predicting a base recirculation flow field with 
shock waves occurring as smeared shocks over several grid points. However, 
such a vast number of grid points might be required to adequately define the 
shock wave patterns that this technique could be wholly unattractive from a 
computational standpoint. It is also questionable if the amount of mass 
reversed by a detached smeared shock near a reflection plane or centerline 
will adequately approximate the mass reversed by a true detached shock which 
will be infinitesimally thin. For these reasons it is recommended that some 
future effort be directed toward improving numerical methods so that shock 
waves could be located between grid points and treated as discrete discontinuities. 
Independent Research efforts at Lockheed/Huntsville into such areas as finite- 
element applications to fluid mechanics could provide an attractive alternative 
technique for the solution of such complex three-dimensional flow fields. 

A literature survey also indicated that there exists a lack of quantitative 
information as t what happens in the region of plume impingement and re- 
sulting shock detachment. Such essential information as the dominating 
parameters controlling the shock detachment height and magnitudes of re- 
versed flow components behind a strong shock are still unknown. This 
problem is akin to blunt body flow, and it is believed that much could be 
learned about dominant terms and controlling parameters as has recently 
been done for blunt bodies by Brainerd and Waldrop (Ref. 10). 
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